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Adopting temperature-dependent interaction parameters in the Lennard-Jones potential, the vapor-liquid phase diagram of 
methane was produced using NVT Gibbs Ensemble Monte Carlo technique. Published second virial coefficient data were used 
to fit a simple two-parameter temperature-dependent model for the interaction parameters. The simulations were carried out in 
the temperature range 120-190 K. The critical density and temperature were evaluated using Ising-scaling model. Using the 
temperature-dependent interaction parameters in the simulation has reduced the root mean square deviation by 94.7% 
compared to the temperature-independent interaction parameters. The evaluated critical temperature was enhanced using 
temperature-dependent interaction parameters, whereas the simulations using temperature-independent interaction parameters 
predict a better critical density value. 
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Introduction 


Methane, the major constituent of natural gas, represents a promising alternative energy source. It is also used to calibrate the density 
transducers which are utilized to directly measure the densities of natural gas transferred in pipelines. It can also be added to some 
refrigerants to enhance their properties (Betaouaf, et al., 2014, Fischer, et al., 1984, Li, et al., 2012, Nie, et al., 2018, Petropoulou, et 
al., 2018, Ungerer, et al., 2007, Uribe-Vargas and Trejo, 2005, Vrabec and Fischer, 1996). Knowledge of the thermodynamic 
properties of methane is necessary to processes of liquefaction, separation, and storage. Therefore, it is required to provide accurate 
thermodynamic properties of this alkane. Among several models available for modeling the intermolecular forces, Lennard-Jones (LJ) 
potential have shown to predict accurate results of certain properties of methane (Fischer, et al., 1984, Jorgensen, et al., 1996, Murad 
and Gubbins, 1978, Saager and Fischer, 1990, Skarmoutsos, et al., 2005, Stassen, 1999); therefore, it has been adopted in this study. 
Saager and Fischer (SF) (Saager and Fischer, 1990) were able to predict quite accurate simple thermodynamic data of liquid methane 
using the interaction parameter values 149.9 K and 3.733 A for the well-depth (€/k) and collision diameter (ø), respectively. These 
values were obtained from fitting vapor pressure of methane and liquid densities. Murad and Gubbins (Murad and Gubbins, 1978) 
have used more complicated five-centered LJ potential to predict more accurate thermodynamic properties of methane. Tchouar et al. 
(Tchouar, et al., 2004) have studied the properties of liquid methane at low temperatures and over a large range of pressure using 
molecular dynamic simulation. They utilized a Feynman-Hibbs temperature-dependent potential form which lead to more accurate 
thermodynamic properties. None of the above-mentioned studies have used temperature-dependent interaction parameters (TDIP) 
model in their simulations. In a previous investigation (Al-Matar, et al., 2008, Al-Matar, et al., 2015), it has been shown that the usage 
of TDIP yields to more accurate vapor-liquid equilibrium properties for Argon. In this study, the vapor-liquid equilibrium of methane 
will be predicted using TDIP and compared with the phase diagram obtained using TIIP, and other modified parameters in the 
literature. 


1 Materials and Methods 


1.1 Intermolecular potential 


The Lennard-Jones intermolecular model is used to represent the pair interaction between methane molecules i and j separated by the 


distance r;j: 
Gij 6 
ace ” 


The second virial coefficient B was utilized to determine the models of the interaction parameters, ¢;; and oj;, as a function of 
temperature. An optimization process was carried out aiming to minimize the least-squares of the errors between the experimental 
values of the second virial coefficient (Dymond and Smith, 1980) and those calculated using the equation: 


oy\ 2 
u(r) = 4Eij (9 
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Betassical = —2mN, f [eo — 1] r?dr (2) 


where N, is the Avogadro’s number and k is the Boltzmann constant. This nonlinear optimization problem was solved using the 
Marquardt-Levenberg (Marquardt, 1963, Press, et al., 1992) method by minimizing the root mean square deviation (RMSD) value, 
given by: 


RMSD = a EG z Bea a 


where M is the number of observations, Bi exp and Bica are the experimental and calculated second virial coefficient for the ith 


observation, respectively. The integral was evaluated using a 21-point Gauss—Kronrod quadrature to estimate the integral and the 
associated errors. 


Only two-parameter relationships between the independent variable, the temperature, and the dependent variables €/k and ø were 
attempted. These relationships involved combinations of linear, reciprocal, exponential and reciprocal—exponential terms. Only 
models with parameters that cause ø to increase and £/k to decrease with temperature were considered feasible (Al-Matar, et al., 2008, 
Al-Matar, et al., 2015). The feasible models were subjected to a model discrimination process based on calculating posterior 
probability (Marquardt, 1963). Table 1 shows the values of TIIP in the literature and the chosen model of TDIP. 


Table 1 Temperature independent values and Temperature dependent 


ee . models of the interaction parameters. 
Within the temperature range of interest for methane 120-190 K, P 


used in the simulations which were carried out, the value of o e/k (K) o (A) 


changes between 3.74 and 3.91 A, and the value of ¢/k changes in 


Temperature Independent Value 


the range 164.75-157.30 K. While the range of ø is in general (Stassen, 1999) 148.2 3.82 
agreement with available temperature-independent values, the values l f 
of €/k are higher by 6.1-11.2%. 

8 y. Temperature dependent Model ones reer 


1.2 Simulation details TCA OTARE oibac GOLY 


Gibbs ensemble is utilized to simulate the coexistence of vapor and liquid phases using systems of 500 atoms. All simulations were 
carried out with a spherical potential truncation for separations greater than 2.5 ø and tail corrections included. The vapor-liquid 
simulations were started using two boxes with simple cubic lattice. The number of atoms was equally distributed between the two 
boxes. The temperature of the entire system is maintained constant and surface effects are avoided by placing each box at the center of 
a periodic array of identical boxes. The simulations were carried out with 200 equilibration and 2000 production cycles. The type of 
Monte Carlo moves were selected at random according to the following probabilities: 0.9089 translation, 0.0909 particle swap and 
0.0002 volume exchange with one volume exchange per cycle. 


The results obtained from the simulations include energy, pressure, density, chemical potential and number of atoms for both vapor 
and liquid phases. The uncertainties in the ensemble averages were calculated by dividing the post-equilibrium results into ten blocks 
then taking the grand average (Frenkel and Smit, 1996). The code employed in this work is developed in-house using object-oriented 
programming in FORTRAN-90 and it is available upon request. 


1.3 Critical constants 


Gibbs ensemble simulations become unstable near the critical region due to the small difference of the free energy between the two 
phases. According to the scaling law, in the case of Gibbs ensemble simulation, the critical temperature is evaluated by the calculated 
p —T coexistence data (Panagiotopoulos, 1994): 


PL — py =b,- T)” (4) 


Where y is the non-classical 3D Ising critical exponent (y = 0.325). The parameter b and the critical temperature T, will be obtained 
from the fit. Subsequently, the critical density p, can be determined using a fit based on the law of rectilinear diameters 
(Panagiotopoulos, 1994): 


L = p, + A(T, — T) 6) 


where the parameter A and the critical density p, will be evaluated from the fit. 
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2. Results and Discussion 


The second virial coefficient data was fitted using TIP as well as TDIP forms. Figure 1 shows the residuals for the second virial 
coefficient data of methane using the values available in the literature o = 3.82 A and e/k = 148.2 K (Stassen, 1999), and fitted 
temperature-dependent models, o = 4.199 — 55.168/T and €/k = 144.54 + 2425.014/T. Using TDIP the model yield to closer 
values of the second virial coefficient than TIP form. The improvement appears in Figure 1 as the difference between experimental 
and calculated values decreases in temperature-dependent case and gets closer to the zero line than the temperature-independent case. 
However, similar to the observation in a previous work (Al-Matar, et al., 2008, Al-Matar, et al., 2015), the residuals still show 
systematic trends suggesting that the temperature-independent parameters are inadequate to describe the behavior of the second virial 
coefficient over a wide temperature range. 


Figure 2 compares the simulation results of the vapor-liquid coexistence curve of this work using a TDIP model and the simulation 
values obtained using TIP taken from the literature. It also shows the simulation results of Skarmoutsos (Skarmoutsos, et al., 2005) 
using Saager-Fischer (Fischer, et al., 1984, Saager and Fischer, 1990) model of LJ potential as well as the experimental data 
(Kleinrahm, ef al., 1988). A summary of the results using the TDIP model is given in Table 2. 


Table 2 Coexistence energies, pressures, and densities, of pure methane for different temperatures at a total number of atoms N=500 


Energy (J/mol) Pressure (bar) Density (g/cm?) 
Temperature (K) Liquid Vapor Liquid Vapor Liquid Vapor 
120 -8331.7 -132.5 9.1 1.34 0.428 2.391x107 
140 -7398.3 -184.8 1.0 3,15 0.372 4.782x107 
160 -6506.4 -404.0 1.7 9.30 0.324 1.339x107 
180 -5612.2 -867.8 21.7 22.0 0.279 3.347x107 
190 -5205.9 -1797.4 40.93 32.44 0.259 6.873x107 
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Fig. 1 Residuals between experimental second virial coefficient and Fig. 2 Vapor-liquid phase diagram curve for methane. (O) Experimental 
calculated second virial coefficient. (O Temperature independent data(Kleinrahm, et al., 1988), (....... ) temperature independent 
case, O Temperature dependent case). curve, (-.-.-.-.-)Skarmoutsos-Saager-Fischer model (Skarmoutsos, et 

al., 2005}, ( ) Temperature dependent curve. ® Experimental 


critical point, W temperature independent critical point, W 
Skarmoutsos-Saager-Fischer model-critical point, ® temperature 
dependent critical point. 


The results of this work obtained using TDIP are closer to the experimental data than the values produced using TIP. The RMSD 
values were calculated to be 5.6 x 107? and 1.07 g/cm? for simulations using TDIP and TIIP, respectively (see Table 3). The gain in 
accuracy as measured by the RMSD is, therefore, about 18 folds compared to the TIIP values in the literature. The deviation between 
the experimental values (Kleinrahm, et al., 1988) of the critical density of methane and the predicted results using TDIP and TIIP in 
the literature are -7.4% and 9.2%, respectively. The vapor-liquid coexistence curve of this work, however, is close to the predicted line 
by Skarmoutsos et al., who adopted the Saager-Fischer (SF) model of LJ potential (Skarmoutsos, et al., 2005). The RMSD obtained 
using the results of Skarmoutsos et al., (Skarmoutsos, et al., 2005) was found to be 4.1 x 10-2 g/ cm?. The RMSD values were 
calculated from the experimental data (Kleinrahm, et al., 1988) of the density in the temperature range 120-190 K. Saager-Fisher 
model of LJ potential was obtained from fitting vapor pressures and liquid densities; and therefore, it is anticipated that the simulations 
of Skarmoutsos (Skarmoutsos, et al., 2005) produce more accurate values in the liquid region. This procedure has also enhanced the 
prediction of the critical density (p,). The values of p, are calculated to be 0.151 g/cm?, with 6.9 % deviation from the experimental 
value, and 0.162 g/ cm, with 0.5 % deviation from the experimental value, using TDIP and TIP of SF model(Skarmoutsos, et al., 
2005), respectively. On the other hand, the predicted critical temperature using TDIP (T, = 193.1 K) was more accurate than reported 
temperature by Skarmoutsos et al. (T, = 194.3 K). 
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Figure 3 illustrates the vapor pressure versus temperature curves obtained from vapor phase results using TDIP and compares it with 
the experimental vapor pressure data, as fitted by the Antoine 
equation, for methane. The values of RMSD in vapor pressure of 
methane using TIIP, TDIP and TIIP-SF model are shown in Table 4. 


Table 3 Comparison between critical densities and temperatures for 
different studies 


The values of RMSD obtained from this work is the highest compared ; F 
to the results using TIIP and SF model. It is anticipated that using the Pe (kg/m*) T(K) RMSD (kg/m*)) 


vapor pressures to fit the LJ potential in the SF model enhances the 
results of vapor pressure predictions. This is due to fitting the 
parameters to the vapor pressure of methane. However, compared toa THP 178.2 190.0 1.07 
previous investigation (Al-Matar, et al., 2008, Al-Matar, et al., 2015), 


Experimental 162.6 190.6 - 


$ . * -2 
the resulting RMSD value in vapor pressure of argon, 6.21 bar is close TIP-SF model 161.8 193.5 41x10 
to the current RMSD result of methane, 6.57 bar. This indicates that TDIP (this work) 151.4 194.3 57x107 
the error in our simulation in vapor pressure has a similar trend in both 
studies. 


Table 4 RMSD in vapor pressure of methane for cases of THP, TIIP-SF 


j ' j j model and TDIP 
50t Temperature dependent 
o Skarmoutsos et al. 

40} Temperature independent | RMSD (bar) 
5 Experimental data THP 0.742 
a 
g 307 ] TIIP-SF model (Skarmoutsos, et al., 2005) 0.830 
= | 
n 
E 207 1 TDIP (this work) 6.568 

10} | 
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Fig. 3 Vapor pressure of methane versus temperature for different systems. 
Conclusions 


Two parameter models of the interaction parameters, 0 = 4.199 — 55.168/T and ¢/k = 144.54 + 2425.014/T in LJ potential were 
used to simulate the vapor-liquid phase diagram of methane. It was shown that using TDIP improves the simulated vapor-liquid 
coexistence curve of methane modeled as an LJ fluid. Adopting TDIP in the simulations reduces the RMSD of density by 94.7% 
relative to the values generated using TIIP. However, the predicted critical values were somewhat mixed; as the calculated critical 
temperature was improved, while the opposite was true for the evaluated critical density. The results of this study are close to the 
simulations of TIP using SF model, which has been obtained by fitting the vapor pressure and liquid densities of methane. Their 
results show closer values to the experimental data notably in the liquid phase region. 


Nomenclature 

A =Ising scaling law parameter [-] 

b =Ising scaling law parameter [-] 
Baassicat =Second virial coefficient [cm?/mol] 
Bica =Calculated second virial coefficient [-] 
Bi exp =Experimental second virial coefficient [-] 

k =Boltzmann Constant [-] 
LJ =Lennard-Jones [-] 
M =Number of values [-] 
N =Number of atoms [-] 
Na =Avogadro’s Numbers [-] 

r =Distance [Å] 
RMSD =Root mean square deviation [-] 
SF =Saager-Fischer [-] 

T =Temperature [K] 
To =Critical temperature [K] 
TDIP =Temperature dependent interaction parameters [-] 
THP =Temperature independent interaction parameters [-] 

u =Interatomic Potential [J] 
fy =Well-depth [cm] 
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y =The non-classical 3D Ising critical exponent [y = 0.325] 
Pe =Critical density [g/cm?! 

PL =Density of liquid phase [g/cm°] 

Pv =Density of vapor phase [g/cm°] 

Gij =Collision Diameter [A] 
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